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The radiation damping effect on the diamagnetic relativistic pulse accelerator (DRPA) is stud- 
ied in two-and-half dimensional Particle-in-Cell (PIC) simulation with magnetized electron-positron 
plasmas. Self-consistently solved radiation damping force converts particle energy to radiation en- 
ergy. The DRPA is still robust with radiation, and the Lorentz factor of the most high energy 
particles reach more than two thousand before they decouple from the electromagnetic pulse. Re- 
sulted emitted power from the pulse front is lower in the radiative case than the estimation from 
the non-radiative case due to the radiation damping. The emitted radiation is strongly linearly 
polarized and peaked within few degrees from the direction of Poynting flux. 



I. INTRODUCTION 

Radiation loss and damping can become important 
in the plasma energetics and dynamics when charged 
particles suffer extreme acceleration. In the ultra- 
relativistic regime, the accumulated effect from ra- 
diation damping can severely limit individual parti- 
cle acceleration even if the radiation damping force 
is weak compared to external forces. However, con- 
ventional Particle-in-Cell simulations of collisionless 
plasmas have not included radiation effects. We de- 
veloped a new 2-1/2-D code including self-consistent 
radiation damping, and studied a effect of radiation 
damping in particle acceleration driven by relativistic 
pulse accelerator (DRPA) in the strong magnetic 
field limit. From the radiation-damped plasma and 
field evolution, we obtained the observable high en- 
ergy radiation output. This radiative PIC simulation 
code is applicable to to a wide range of high-energy 
astrophysics phenomena (pulsars, blazars, gamma-ray 
bursts) and ultra-intense laser applications 0- 

The outline of this paper is as follows. In Sec. [H] 
the derivation of the radiation damping force in the 
relativistic form is shown. In Sec. IIIII we explain how 
the radiation damping force is implemented into the 
PIC code. Results of simulation and comparison with 
non-radiative case are given in Sec. II VI and we sum- 
marize in Sec. 



II. DERIVATION OF RADIATIVE FORCE 

In the conventional PIC simulations, radiation 
damping force f rad is ignored because of its small am- 
plitude compared to external forces F ext . However, it 
is not negligible when the accumulated work done by 
the radiation damping during deceleration time Tdecl 
is comparable to the work done by the external force 
in typical acceleration time . 



It is impractical to include high-frequency radiation 
wave into the electromagnetic calculation in the PIC 
code, since the typical radiation wavelength is much 
shorter than the spatial resolution of the fields (~ De- 
bye length Xd = c/uj pe , where u> pe — Ait pe / m e is the 
electron plasma frequency) . Accelerated particles can 
emit up to the critic al frequenc y uj c = 3j 2 fl ce , where 
7 = E/rn e c 2 = l/^/l — v 2 / c 2 is the Lorentz factor 
and Q C e = eB/(m e c) is the electron gyro-frequency. 
The ratio of the critical radiation wavelength A c to 
Xd is given by X c /Xd = {2Truj pe )/ (3~f 3 fl ce ), which is 
<C 1 because w pe /f2 ce < 0.1 in magnetic-dominated 
cases and 7>1. 

Instead, we introduce a radiation damping force in 
the form of the Dirac-Lorentz equation 3] , which is 
proven to be exact for a classical point particle Q. 
The 4- vector form of Dirac-Lorentz equation with the 
damping force g % is given by 
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F ik u k +g l 
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where F is the electromagnetic tensor ; 5], u % is the 
velocity four-vector, and g l is the radiation damping 
force 
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which satisfies the auxiliary relation g % Ui = 0. Since 
Eq. (J2J) includes the second derivative of u 1 , there ex- 
ist unphysical run- away solutions 5j. To avoid the 
run-away solutions, we assume gi <C eF lk Uk/c and 
eliminate the second derivative terms in g % by express- 
ing d 2 u % jds 2 with the first derivative of Eq. with 

a 1 = o, 
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Using the fact (dF lk /dx l )uiUk = 0, we find 
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The resulting damping force term f rad in the 3- vector 
form is given by [3j 

frad — 777^ krad X 
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where d is the 3-velocity, and E and B are the self- 
consistent electric and magnetic fields. Here we intro- 
duce a non-dimensional factor k rad given by 

k rad = = 1.64 x 1CT 16 x ^(gauss), (6) 

c 

where r e — e 2 /(mc 2 ) is the classical electron radius. 

The first term of the radiation damping force ijHJ) 
represents the radiation damping due to the pondero- 
motive force acceleration. The third term is Compton 
scattering by large scale (A > A_d) electromagnetic 
field which reduces to Thomson scattering in the clas- 
sical limit 6]. Note here that the scattering between 
radiation field and particles is not considered since 
electromagnetic field in Eq. © is averaged over the 
grid separation that is much larger than the wave- 
length of radiation, as we discussed above. In other 
words, plasma is perfectly transparent to emitted ra- 
diation, and all the radiation disappears from the sim- 
ulation box immediately after particles radiate. 

The magnetic field strength of pulsars are in the 
range of 1CT 4 < k rad < 10~ 3 , and k rad ~ 10~ 2 
for magnetars. In the laboratory, k rad ~ 10 -6 
(I/10 22 Wcm -2 ) 1 / 2 for lasers of intensity /. Hence 
for ultraintense laser interactions (/ > 10 22 Wcm~ 2 ), 
radiation damping effects can become significant on 
simulation time scales > 0.3ps (= lO 5 ^^ 1 ). (See 
also 2J) In our simulation, we can enhance the radia- 
tion effect by increasing the value of k rad , or the 'effec- 
tive' electron radius. However, we should restrict our- 
selves not to reach the quantum- limit, hfl ce ~ m e c 2 
or B > 4.4 x 10 13 gauss, which corresponds to k rad = 
7.2 x 10~ 3 , or the Dirac-Lorentz equation Q fails. We 
choose k r ad from zero to 10 -3 in the simulation to en- 
hance the radiation effect and \f ra d\ T sim — \F e xt\^ce 
so we can see the difference between radiative and 
non-radiative (k rad = 0) case within the simulation 
time-scale T S i m ~ lO 4 ^" 1 . 



III. IMPLEMENTATION OF RADIATION 
FORCE 

The 2-1/2D explicit PIC simulation scheme is used 
with the explicit leap-frogging method for time ad- 
vancing 0. Spacial grids for the fields are uniform in 



both x and z directions, Ax = Az = Xo- The simu- 
lation domain in the x— z plane is —L x /2 < x < L x /2 
and < z < L z with a doubly periodic boundary 
condition in both directions. 

Following Liang et. al. 0, the initial plasma is uni- 
formly distributed at the center of the simulation box, 
— 6Aa; < x < 6 Ax and < z < L z . The background 
uniform magnetic field Bq = (0, Bq, 0) is applied only 
in the same region, so that the magnetic field freely 
expands toward the vacuum regions, x > 6Ax and 
x < —6Ax with accelerating plasma. We choose L x 
to be long enough so that plasma and EM wave never 
hit the boundaries in the x direction within the simu- 
lation time. 

The initial temperature of plasma is assumed to be 
a spatially uniform relativistic Maxwellian, ksT e = 
ksTp — lMeV, where the subscripts e and p refer to 
electrons and positrons. 

The radiation damping force JSJ is calculated self- 
consistently and fully-explicitly as follows. The veloc- 
ity of each particle should be updated each time step 
from v(t — At/2) to v(t+ At/2), using the electromag- 
netic field and f rad at time t following the equation 
of motion in the relativistic form, 



dp 
~d~t 



eE+-vxB + /,,, 



(7) 



where p = mjv is the relativistic momentum. All the 
terms in the equation JSJ are given by the ordinary 
leap-frog field solver except dE(t)/dt term, which is 
not calculated at time t in the conventional PIC code. 

In order to calculate dE(t)/dt, we need to know 
J(t). First, we advance the position x half timestep 
from t — At/ 2 to t for each particle using velocity 
v(t — At/2). Next, temporal current Jt{t) is calcu- 
lated from v(t — At/2) and x(t). Note that this cur- 
rent J t (t) is not exact because of the velocity field not 
at time t but v(t - At/2). Finally, the term dE(t)/dt 
is calculated using the Maxwell equation, and the ra- 
diation damping force is calculated for each particle. 
To update the velocity, we consider a sum of the radi- 
ation force and the electric field force as a net 'accel- 
eration' force, and apply the Boris method for particle 
gyration 0- 



IV. RESULTS 

We choose six different sets of parameters shown 
in Table [I] by changing k rad = 0, 10~ 4 , 10~ 3 and 



0.1, 0.01, and run simulations for each case. 
Hereafter, we call Run A-C as the weak magnetic field 
case and Run D-F as the strong magnetic field case, 
based on the ratio uo pe /VL ce . 

First, we check the total energy conservation for 
Run D,E and F in Fig.Q] In the radiative (RD) cases, 
the energy loss by the radiation E rad is obtained from 
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TABLE I: Parameters for each runs 









Duration tQ ce 


Run A 





0.1 


10000 


Run B 


10" 4 


0.1 


10000 


Run C 


1(T 3 


0.1 


10000 


Run D 





0.01 


70000 


Run E 


HT 4 


0.01 


70000 


Run F 


1(T 3 


0.01 


70000 




FIG. 1: System-integrated energy in the electromagnetic 
field (1), particles (2), radiation damping (3), sum of field 
and particle energy (4), and total energy (5) functions of 
time for Run D (solid), E (dashed) and F (dash-dot). 
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the time integral 

Era d (t) = J* V W ■ fraAt')^ dt>. (8) 

In the non-radiative (NRD, k = 0) case, total energy 
of the system is given by the sum of the held energy 
Efi e (Line 1) and kinetic energy E^m (Line 2), and 
it conserves. In the RD case, however, sum of E^m 
and Ef ie does not conserve (Line 4), but sum of Eki n: 
Efi e and the radiation energy E rac i (Line 3) conserves 
(Line 5), indicating that the radiation damping force 
is self-consistently calculated. In all the RD cases, the 
energy is transferred from field to particle, and then 
radiation, indicating the DRPA mechanism acceler- 
ates particles even in the RD cases. Energy transfer, 
however, from field to particles becomes less efficient 
with larger radiation damping force. Radiation pre- 
vents energetic particles to get accelerated from EM 
field, resulting less efficiency of transfer in high energy 
tail of the particle distribution. 

Figure |2 shows the momentum distribution of par- 
ticles for (a) Run A(x < 0) and Run C(a; > 0), and 
(b) Run D(x < 0) and Run F(x > 0). The DRPA ac- 
celerates electrons and positrons in the same direction 



FIG. 2: Momentum distribution of particles for Run A 
[(a), x < 0] and Run C [(a), x > 0] at tQ. ce = 5000 (green, 
blue) and tQ ce = 10000 (magenta, red), and Run .D [(b), 
x < 0] and Run F [(b), x > 0] at tQ ce = 35000 (green, 
blue) and tQ. ce = 70000 (magenta, red). Results in the 
positive and negative x directions are identical in all cases. 



along the x axis, whereas electric field E z accelerates 
electrons in the positive z and positrons in the nega- 
tive z direction respectively in the positive x region. 
In the negative x region, acceleration directions are 
opposite for both species, forming X shape distribu- 
tion in the p x — p z plane as a result. Charge sep- 
aration does not occur because of the periodicity in 
the z direction. Resulted induced current J z acceler- 
ates particles in the x direction by the ponderomotive 
force J x B. The ponderomotive force creates succes- 
sive 'potential wells' in the x direction, which captures 
and accelerates co-moving particles. We emphasize 
here that there is no charge separation in the x direc- 
tion because of no mass difference between electron 
and positron. 

Obviously, particle momenta in both x and z di- 
rections are radiated away in the RD cases in both 
weak and strong magnetic field cases. Bifurcation in 
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FIG. 3: Average instantaneous radiation power from par- 
ticles within 30Ao from pulse front for (a) uj pe /Q. ce = 0.1 



and (b) 0.01 with k ra d 



(dash-dot line) and k ra d 



10 4 (dashed line). Estimated power for NRD cases with 
krad = 10" 3 (solid line) and 10~ 4 (dotted line) arc also 
shown for comparison. 



high energy tails occurs in Run C since slow parti- 
cles can not keep up with the speed of the first (and 
the fastest) ponderomotive potential well and decou- 
pled from it. Then the second potential well captures 
these slow particles and accelerates them again, which 
creates the bifurcated high energy tails in the phase 
space. The bifurcation is less clear in Run A than Run 
C because energetic particles are not suffered from the 
radiation damping and most of them can stay in the 
first potential well. 

In the strong magnetic field case [Fig. |2tb)] , accel- 
eration is strongly reduced in the RD case (Run F). 
For high energy (7 3> 1) particles, the Compton scat- 
tering [the third term in Eq. becomes a dominant 
term because of the large Lorentz factor, and makes 
the DRPA acceleration less efficient. Bifurcation is 
strongly suppressed in both Run D and F, since the 
first potential well is deep enough to capture almost 
all of the energetic particles. 

Next, we compare the radiation power of RD cases 
with NRD cases. For the NRD case, we estimate 
the radiation power using the relativistic dipole for- 
mula 01 



(9) 



where F\\ and F±_ are the parallel and perpendicular 
components of the force with respect to the particle's 
velocity. The bracket () indicates that we take the 
average of all the particles locating within 30Xd from 
the pulse front, and the number of particles within this 
layer decreases with time because of the decoupling 
from the EM pulse. Since Eq. (O is proportional to 
krad, we plot Run A and D with k = 10~ 3 and 10 -4 
to compare with the other four RD cases. 

For the RD cases, the radiation power is calculated 
using the formula 



(P) = \frad ■ 



(10) 



We compared the result of these two formulae 
for the RD cases, and they matched within the 
lincwidth. Therefore, We only show the result of Eq. 
ltTU|) here. 

Figure [3] shows the instantaneous radiation power 
(P) for all runs. The initial peak around tQ, ce = 100 
is due to thermal-cyclotron emission |6(], and the es- 
timated radiation power for the NRD case quantita- 
tively matches with the RD case. At a later time 
(tQ ce > 1000), however, more power is irradiated in 
the NRD cases than in the RD cases, because energetic 
particles are continuously accelerated without losing 
their energy self-consistently in the NRD cases. In the 
RD cases, however, all the energetic particles are de- 
celerated by the radiation damping, and the instant 
radiation power decreases caused by slower velocity 
and smaller damping force. 

Finally, we calculate the radiation field and its an- 
gular dependence self-consistently using the velocity 
and acceleration of each particle. Intensity / and po- 
larization II of the radiation received by the observer 
located at x are given by p|, |(| 



ret 



(11) 



and 



E 2 y (n,r) 



Ei ■ y\ 2 



ret 



ret 



U(h,r) = 2j2m-y)(Ei-z)] Tet , (12) 

i 

{Elf + (£2)2 - 2E%E* + C/ 2 



n(n,r) 
where 



m + E i 



Ei = 



e nx [(n-Pj) x_/3j 



-,(13) 



(14) 



is the radiated electric field from particle i located at 
r, n is a unit vector in the direction of x — r(r), (3 = 
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FIG. 4: Contour plots of instantaneous intensity log 10 I 
as a function of local time t for 9 — cj> — for (a) Run C 
and (b) Run F, and intensity (red lines, right scales) and 
polarization (blue lines, left scales) as functions of obser- 
vational time t for (c) Run C and (d) Run F. Intensity is 
in arbitrary scale. To obtain time dependence of intensity, 
instantaneous intensity from each particle is summed up 
along the light cone r = t— R/c =const., shown as a black 
dotted line in panel (a). The light cone moves horizontally 
leftward with t. 



v(t)/c, and /3 = df3/dt. We assume that \x\ ^S> \r\ so 
that h is parallel to x. The square brackets with a 
subscript "ret" mean that the quantity in the brackets 
is evaluated at the retarded time r = t — R/c, where 
R = |x — r\. To specify the direction of the observer 
with respect to the x axis, we introduce and <p as 



(cos 9 cos cj), cos sin (f>, sin 9). 



(15) 



In Figs.Qfa) and (b), the contour plot of instanta- 
neous intensity before taking the summation over the 
retarded time is plotted, as a function of local time 
t with <fi = 9 = 0, illustrating the ray-tracing tech- 
nique used in Eqs. (|llfl and (|12f> . We take a sum of 
intensity along the light cone r = t — R/c =const., 
which is indicated as a black line, up to tVl ce = 10000 
for Run C and m ce = 70000 for Run D. The light 
cone moves toward the negative x direction with r, 
and we take r = when the pulse front reaches to the 
observer. Since all the particles cross the light cone 
r =const. only once, we take a sum of / and II as fol- 
lows. First, we follow the trajectory of each particle, 
and store their instantaneous intensity. When each 
particle crosses the light cone rf2 ce = 1, we take the 
average of the instantaneous intensity of the particle 
over the trajectory between rfl ce = and 1, and add it 
to the intensity /. Then we continue to follow the par- 
ticle trajectory and sum up the intensity again until it 
hits the next light cone. We continue the summation 
until the signal in the initial plasma thickness reaches 
to the observer, r = 12A rj/c, We also note that the 
intensity is extremely asymmetric because only ener- 
getic particles accelerated in positive x direction can 
radiate strong emission to the observer. 

The time dependence of intensity and polarization 
with cf> = 9 = is shown in Figs. 0fc) and (d). If 
all the particle moves with the speed of light and are 
continuously accelerated toward the direction of n, 
all the radiation emitted from particles should reach 
the observer at the same time and the signal becomes 
a <5-function pulse. However, energetic particles are 
bouncing back and forth within the ponderomotive 
potential well, and slower particles are dropped off to 
the next well, which broaden the spatial distribution 
and the resulting radiation duration. 

Intensity / is shown as solid lines in Fig. 0Jc) and 
(d), indicating the duration time is 20rf2 ce for Run 
C and 200 ~ 400rS! ce for Run D with single peak. 
Since the initially applied electromagnetic field is lin- 
early polarized, we expect that the radiation is also 
strongly polarized, with the small depolarization com- 
ing from the initial random velocity distribution in the 
y direction. Polarization II is shown as dotted lines 
in Fig.^Jc) and (d), illustrating that the radiation is 
strongly linear-polarized as anticipated. 

In Fig. [SI we show the total radiation intensity 
J I(cf>,9)dT. The radiation is a very short pulse in 
Fig. ^because we consider single simulation box cen- 
tered at the origin only. However, because of the peri- 
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FIG. 5: The total radiation intensity J I(4>,8)dT for Run 
C and Run F as a function of the angle (f>, 9). 



odic boundary condition in the z direction, we should 
consider the multiple simulation boxes along the z 
axis, and consider time delay and angle difference to- 
ward the observer from each simulation box, which is 
very complex even in the Cartesian coordinate. 

Instead of including these geometrical effects into 
the intensity calculation, we simply compare the total 
intensity integrated over t as a function of the angle. 
Intensity peaks around 6> = 3^8° in = cases 
(solid and dash-dot lines), corresponding to the di- 
rection of high energy particles in Fig. [5] Intensity 
rapidly decreases with both <p and 9, indicating radi- 
ation is strongly collimated in the x direction. 

Intensity distribution in the z direction is due to 
the initial electric field acceleration, as we discussed. 
Since there is no acceleration in the z direction after 
the induced current is formed, eventually all the mo- 
mentum in the z direction will be emitted away, which 
narrows the intensity distribution toward positive x 
direction. At the same time, however, decoupling of 
particles from EM pulse in the x direction makes the 
DRPA less efficient in the x direction, which widen the 



distribution. Thus, the intensity in the x — z plane is 
always distributed over a finite angle, but the distri- 
bution in much later time is still an open question. 



V. SUMMARY 

In summary, we observed the self-consistent radia- 
tion damping effect on the interaction of the DRPA 
with electron-positron plasma via a relativistic PIC 
simulation. We have found that field and particle 
energies are transferred to radiation, and the cou- 
pling between the field and particles becomes less effi- 
cient with larger radiation damping. Comparison with 
the non-radiative case showed that radiation damping 
force decelerates the energetic particles accelerated by 
the DRPA, and resulting radiation power is smaller 
in the radiative case. The radiation field is strongly 
linearly polarized both in weak and strong magnetic 
field cases, which may be detectable by 7-ray burst 
observations or laser experiments as an indication of 
the DRPA mechanism. The simulations shown here 
are too short to see multiple peaks in an intensity 
time profile as seen in GRB observations, and to de- 
termine the final radiation pattern after particles are 
completely decoupled from EM pulse. These two ques- 
tions remain as future problems. 
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